Logistic Regression
James Rigby
10/26/2019
In this tutorial, I am going to describe the inner working of logistic regression. The first portion of the post is going to break the logistic model into small elements, writing a function for each of the chunks. This is meant to illustrate how logistic regression is actually estimated. The second part is going to approach logistic regression from a bayesian perspective and illustrate how Bayesian learning can be implemented to create an online machine learning tool. The final portion of this paper is going to demonstrate how to speed up bayesian analyses using approximations to the posterior. This post is going to be heavy on statistical theory, but I will try to point out the sections that those primarly conserned with application should be aware of in the section headers.
Problem Setup
Your consulting firm was recently hired to build a system that integrates with the company’s Buisiness Intelligence (BI) software to understand the factors that influence employee recruitment. The primary critereon of interest is job acceptances - candidates decisions to accept or reject an offer made my the organization. You decide to begin this project by exploring a logit model (logistic regression), because you hear that it is easily scaleable and able to be developed into an active learning.
The BI platform provides you with information from recent job analyses, salary market surveys, job market research, and the benefits packages that were offered to the employee. I print the head of the data frame below.
Simulating Data
| job_id | applicant_id | salary | job_alt | difficulty | flexibility | market_sal | decision |
|---|---|---|---|---|---|---|---|
| 12 | 1 | 25.00807 | 1 | 1 | 1 | 23.42748 | 1 |
| 39 | 2 | 28.52853 | 2 | 2 | 1 | 27.61328 | 1 |
| 8 | 3 | 27.26684 | 1 | 4 | 1 | 21.40710 | 0 |
| 37 | 4 | 19.09980 | 4 | 3 | 4 | 17.79660 | 1 |
| 1 | 5 | 23.96467 | 1 | 2 | 1 | 24.36497 | 0 |
| 49 | 6 | 27.38086 | 4 | 4 | 3 | 18.26968 | 1 |
Approaching Classification from a Linear Perspective
Why is linear regression not the most appropriate technique for categorical outcomes? If we code job acceptance as 0 or 1 we can predict it using a linear regression model. Furthermore, if the goal of our analysis make predictions about those who accept a job and those who do not accept a job we could create a decision boundary. For example, we could classify job applicants with a predicted value above .5 as accepters and classify those with a score below .5 and rejecters.
p<-dec%>%
ggplot(aes(x = salary, y = decision))+
geom_smooth(method = "lm", se = FALSE, color = "red", fullrange = TRUE)+
geom_point()+
geom_vline(xintercept = 25.80, lty = 3)+
labs(title = "Linear Approach to Classification",
x = "Salary",
y = "Predicted Values (1 = Accepted Job Offer)")+
scale_y_continuous(breaks = seq(0, 1, by = .1))
plotly::ggplotly(p)p<-p+
geom_point(data = data.frame(salary = 55, decision = 1), aes(x = salary, y = decision), shape = 3)+
scale_y_continuous(breaks = seq(0, 1.2, by = .1))## Scale for 'y' is already present. Adding another scale for 'y', which
## will replace the existing scale.plotly::ggplotly(p)Logisic Regression: A Probablistic Approach to Classification
We can still approach classification from a linear perspective, however, we need to transform the function to accomidate the binary data. Logistic regression does this by assuming the log odds of the categorial outcome is linearly related to the independent variables.
\[log\frac {p(y=1)}{p(y=0)} = X\beta+e\]
We can rearrange this function to depict the model in probablistic terms.
\[\frac {p(y=1)}{p(y=0)}= e^{X\beta+e}\] \[\frac {p(y=1)}{1-p(y=1)}= e^{X\beta+e}\] \[p(y=1)= (e^{X\beta+e}-e^{X\beta+e}p(y=1)) \] \[p(y=1)+e^{X\beta+e}p(y=1)= (e^{X\beta+e}) \] \[p(y=1)(1+e^{X\beta+e})= (e^{X\beta+e}) \] \[p(y=1) = \frac{(e^{X\beta+e})}{(1+e^{X\beta+e})} \] This can be written in a slightly more compact form
\[p(y=1) = \frac{1}{(1+e^{-(X\beta+e)})}\] Notice that there is still a linear component \(X\beta + e\), however it is transformed to ensure that the function is bounded by 0 and 1.
How does this new function look? Unlike the linear classifier (red line) the logistic model’s (blue line) predicted values can be directly interpreted as the probability of an employee accepting a job given the predictors.
p<-p+
geom_smooth(method = "glm",
method.args = list(family = "binomial"),
fullrange = TRUE,
se = FALSE,
color = "blue")+
scale_y_continuous(breaks = seq(0,
1.2,
by = .1))+
labs(title = "Comparing Linear and Logistic Approaches to Classification")## Scale for 'y' is already present. Adding another scale for 'y', which
## will replace the existing scale.plotly::ggplotly(p)Zooming out provides a better illustration of how the two functions differ. The logit model assymptotes at 0 and 1 while the linear model is, well, linear. Interstingly, the decision boundaries for the two models seem aligns at a .5 threshold. However, if we adjusted the decion boundary above or below .5, the two models would diverge.
p<-p+
lims(x = c(-20, 70))+
geom_hline(yintercept = 1, lty = 3)+
geom_hline(yintercept = 0, lty = 3)
p$layers[[3]]<-NULL
plotly::ggplotly(p)Lets define the functional form of the logistic regression. To make it’s relation to the linear model explicit, I defined a function called v which is the estimate of the exponentiated portion of the model (i.e., \(X\beta\)). I then define a function called p_hat which transforms the linear output to estimate probabilities.
# Linear combination of predictors
v<-function(X, beta){
v<-X%*%beta
v
}
# Estimated Probability through sigmoid transformation
p_hat<-function(v){
p<-1/(1+exp(-v))
p
}Maximum Likelihood Estimation
Just like linear regression selects the regression coefficients that minimize the sum of squared error, logistic regression has an objective criterea, or objective function. Logistic regression is fit using maximum likelihood estimation. Broadly speaking, maximum likelihood tries to find the parameters (\(\beta\)) that maximize the probability of observing the data that was observed (i.e, \(p(y|\beta, X)\)). The likelihood function for logistic regression is written below.
\[L(\beta) = -\sum_{i=1}^n\ y_i\ ln(\ \hat p_i(y = 1)\ )+(1-y_i)\ ln(\ 1-\hat p_i(y = 1)\ )\]
What does the likelihood function mean? Let’s break the function into two pieces by spliting it at the plus sign.
Part 1
\[y_i\ ln(\ \hat p_i(y = 1)\ )\] This portion of the likelihood function multiplies the oberved outcome (\(y_i\)) by the log of the expected probability. As probabilities fall between 0 and 1 and the log of numbers closer to 0 become more negative, this portion of the likelihood is largest when 1) the predicted probability is high and 2) the observed outcome was one. In other words, this portion of the equation is the largest when observations are correctly classified as belonging to class 1. Observations that are 0 do not contribute to this part of the likelihood. The plot below illustrates how this function behaves across classes as the predicted probability increases.
expand.grid(p_hat = seq(.001, 1, .001), y_i = c(0,1))%>%
mutate(likelihood = y_i*log(p_hat),
y_i = as.factor(y_i))%>%
ggplot(aes(x = p_hat, y = likelihood, color = y_i))+
geom_line()+
labs(x = "Predicted Probability",
y = "Part 1 Likelihood",
title = "Illustration of Part 1 of Likelihood Function")+
guides(color = guide_legend("Observed Y"))+
scale_color_manual(values = c("red", "blue"))
Part 2
\[(1-y_i)\ ln(\ 1-\hat p_i(y = 1)\ )\]
The second half of the likelihood formula follows the same logic. The major differences is that it focuses on encorporating the other class. This portion of the likelihood is maximized when 1) the predicted probability is closest to 0 and 2) the observed class was 0. Thus, this portion is maximized when the model accurately predicts class 0. The plot below depicts this function.
expand.grid(p_hat = seq(0, .999, .001), y_i = c(0,1))%>%
mutate(likelihood = (1-y_i)*(log(1-p_hat)),
y_i = as.factor(y_i))%>%
ggplot(aes(x = p_hat, y = likelihood, color = y_i))+
geom_line()+
labs(x = "Predicted Probability",
y = "Part 1 Likelihood",
title = "Illustration of Part 1 of Likelihood Function")+
guides(color = guide_legend("Observed Y"))+
scale_color_manual(values = c("red", "blue"))
Taken together, the likelihood function is the most positive for the \(\beta\) that do a the best job discriminating between the classes. The function below writes the likelihood function. Its sole arguments are the observed values y and the predicted probabilities p which can be generated by p_hat().
# Implementation of the likelihood
ll<-function(y, p){
-sum(y*log(p)+(1-y)*log(1-p))
}Newtonian Optimization to Solve for Beta
How do we solve for the \(\beta\) that maximize the likelihood of our observed data? One method is to use Newtonian optimization to solve for the maximum of our likelihood function. Newtonian methods rely on a second order Taylor expansion that approximates the likelihood function. So long as the function is second differentiable we can iteratively update our estimates for \(\beta\) using the following equation:
\[\beta_{new}= \beta_{old}-H_\beta^{-1}\Delta\beta\]
where \(\Delta\beta\) is the vector of first derivatives of the likelihood function with respect to \(\beta\) and \(H_\beta\) is equal to the the hessian, or matrix of second derivatives, of the likelihood function with respect to the regression coefficients.
This function sets the derivatives of the quadratic approximation to 0, therby solving for the approximation’s minimum or maximum.
To visually illustrate what newtonian optimization does, Let’s try to optimize the following equation.
\(y = .25x^2\)$
data.frame(x = -20:20)%>%
mutate(y = .25*x^2)%>%
ggplot(aes(x = x, y = y))+
geom_line()
To begin optimizing, we need to initalize the function by picking starting values. Lets pick a point thats really far off from the stationary point (the minimum).
obj_function<-function(x){
.25*x^2
}
data.frame(x = -30:30)%>%
mutate(y = obj_function(x))%>%
ggplot(aes(x = x, y = y))+
geom_line()+
geom_point(aes(x = -20, y = obj_function(-20)))
The derivative of the objective function in this example is equal to \(.5x\) while the second derivative is equal to \(.5\). To verify my calculations are correct, I plot the objective below. As a reminder the first derivative is equal to the tangent line, or the instintaneous slope at a point \(x_n\).
obj_function<-function(x){
.25*x^2
}
data.frame(x = -30:30)%>%
mutate(y = obj_function(x))%>%
ggplot(aes(x = x, y = y))+
geom_line()+
geom_segment(x = -20-3, xend = -20+3, y = obj_function(-20)+10*3, yend = obj_function(-20)-10*3, color = "red", size = 1.07)+
geom_point(aes(x = -20, y = obj_function(-20)))
Let’s take a look at the taylor expansion of this equation about our starting point.
\[f_T(x)= f(x_n)+f'(x_n) \Delta x+ \frac{1}{2} f''(x_n) \Delta x^2\]
We can plug in our starting values and the derivatives to define and plot the quadratic approximation.
\[f_T(-20) = 100 + (-10 \times (x+20)) + \frac{1}{2} (.5 \times (x+20)^2)\]
Plotting this function over the objective function reveals that it is a fantastic approximation to the objective function, which may not be supprising given that true form of the objective is quadratic.
taylor_exp_neg_20<-function(x){
obj_function(-20)+ (x+20)* -10+ 1/2* .5* (x+20)^2
}
data.frame(x = -30:30)%>%
mutate(y = obj_function(x),
taylor = taylor_exp_neg_20(x))%>%
ggplot(aes(x = x, y = y))+
geom_line()+
geom_line(aes(x = x, y = taylor), color = "red", lty = 2, size = 1.07)+
geom_point(aes(x = -20, y = obj_function(-20)))
We can apply the update rule by applying the following formula
\[x_{new} = x_{old}-\frac{\Delta x}{\Delta \Delta x}\] \[x_{new} = -20-\frac{-10}{.5}\]
Given the perfect congruence between the taylor expansion and the objective function, this optimization problem converges in one iteration. This only occurs when the function you are trying to maximize or minimize is perfectly represented by the taylor expansion. Unfortunately, the likelihood function requires multiple taylor expansions to converge but this isn’t a huge deal. We just need to know a little bit of programming to get the job done!
data.frame(x = -30:30)%>%
mutate(y = obj_function(x),
taylor = taylor_exp_neg_20(x))%>%
ggplot(aes(x = x, y = y))+
geom_line()+
geom_line(aes(x = x, y = taylor), color = "red", lty = 2, size = 1.07)+
geom_point(aes(x = -20, y = obj_function(-20)))+
geom_segment(aes(x = -20, xend = -20, y = obj_function(-20), yend = obj_function(-20+10/.5)), lty = 2, color = "red")+
geom_segment(aes(x = -20, xend = -20+10/.5, y = 0, yend = 0), color = "red", lty = 2)+
geom_point(aes(x = 0, y = 0), shape = 2)
Because the likelihood in logistic regression is second differentable, we can apply newtonian methods! The first derivative of likelihood function can be estimated using the following formula
\[\Delta\beta = X^\top (y-p(y=1))\] I implement this function below.
deriv<-function(X, y, p){
t(X)%*%(y-p)
}The second derivative has quite a special place in psychological and assessment methodology. It is referred to as the information matrix. It’s scalar representation may seem a little bit tricky, but the matrix calculations (how I implement my R code below) is actually quite easy!
\[I(\hat\beta) = \frac{\delta^2L(\beta)}{\delta\beta_j\delta\beta_l} = -\sum_{i = 1}^n\ x_{ij}\ x_{il}\ p_i(y = 1)\ (1-p_i(y = 1))\] This can be solved using the following matrix algebra.
\[I(\hat\beta) = \frac{\delta^2L(\beta)}{\delta\beta_j\delta\beta_l} = -X^\top S X\] where S is a diagnol matrix with \(p(y=1)(1-p(y=1))\) along the diagnol.
# Diaganol Matrix of S
S_mat<-function(p){
var<-p*(1-p)
dim<-length(var)
s<-matrix(0, nrow = dim, ncol = dim)
diag(s)<-var
s
}
# Calculate Information Matrix
information_mat<-function(X, S){
-t(X)%*%S%*%X
}Iteratively Reweighted Least Squares
If you read much literature on logistic regression it is only a matter of time before you come across the term “Iteratively Reweighted Least Squares”. An interesting property of the Newton step described above is that it simplifies to another optimization method called iteratively reweighted least squares. This method is explained elsewhere but for a fantastic and detailed description linking the two optimization methods see the stacks exchange post by jld here.
Formulaically, the iterative least squares step is equal to
\[I(\hat \beta)^{-1}X^\top S z\]
where z is defined as \(X \beta+S^{-1}(y-p(y=1))\).
# Define z
z<-function(X, theta, S, y, p){
X%*%theta+solve(S)%*%(y-p)
}# Update the regression weights
update_beta<-function(i_mat, X, S, z){
(solve(-i_mat))%*%t(X)%*%S%*%z
}When do we stop?
Iterative methods require some stopping criterea. The default is for many software is the relative gradient convergence critereon. I also include a maximum iterations option as a safe gaurd, but typically models converge after a relatively small number of iterations.
relative_grad_crit<-function(deriv, i_mat, ll, tol){
abs((t(deriv)%*%solve(i_mat)%*%deriv)/(abs(ll)+1e-6))
}Putting it all together
Below I create a function that pulls together all of the pieces
log_reg<-function(X, y, max_iter = 20, tol = 1e-8){
# Initalize the number of iterations so max_iter will serve as a stopping critereon
i<-1
# Provide starting values for the regression weights
theta<-rep(0, times = ncol(X))
# Solve all components for starting values
# These just implement our functions above
# Solves for the predicted probability given starting values
v_i<-v(X = X, beta = theta)
p<-p_hat(v = v_i)
# Solves for the first and second derivaties of likelihood function given start values
deriv_i<-deriv(X, y, p)
S<-S_mat(p = p)
i_mat<-information_mat(X = X, S = S)
# Initializes the stopping critereon
grad_crit<-relative_grad_crit(deriv_i, i_mat, ll(y, p))
# Implementing Iteratively Reweighted Least Squares
# Runs until stopping critereon is reached
while(i<max_iter&grad_crit>tol){
# Update beta using IRLS step
z_i<-z(X = X, theta = theta, S = S, y = y, p = p)
theta<-update_beta(i_mat = i_mat, X = X, S = S, z = z_i)
# Solve for VI and P with new values
v_i<-v(X = X, beta = theta)
p<-p_hat(v = v_i)
# Sovles for second derivatives
S<-S_mat(p = p)
i_mat<-information_mat(X = X, S = S)
deriv_i<-deriv(X, y, p)
# update ll
ll_new<-ll(y = y, p = p)
# Return a message with the iteration value
message("Iteration: ", i, ll_new)
# move one iteration forward
i<-i+1
# Calculate stoping critereon
grad_crit<-relative_grad_crit(deriv_i, i_mat, ll(y, p))
}
# Calculate final log likelihood and betas
final_ll<-ll(y = y, p = p)
final_beta<- theta
# Format the betas
names(final_beta)<-c(colnames(X))
# Return the Values
list(betas = final_beta,
LL = final_ll,
vcov = -solve(i_mat))
}Giving the Model A Try
# Prep Matrices
X<-dec[3:7]
y<-dec["decision"]
X<-cbind(1, X)
X<-as.matrix(X)
y<-as.matrix(y)
colnames(X)[1]<-"(Intercept)"
our_reg<-log_reg(X = X, y = y, max_iter = 100)## Iteration: 1121.86552074049## Iteration: 2116.080315519869## Iteration: 3115.545280251892## Iteration: 4115.538739109569## Iteration: 5115.538737957367our_reg$betas## decision
## (Intercept) -0.1935129
## salary 0.3384038
## job_alt 1.0078720
## difficulty -1.8396173
## flexibility 0.2497007
## market_sal -0.2805268
## attr(,"names")
## [1] "(Intercept)" "salary" "job_alt" "difficulty" "flexibility"
## [6] "market_sal"our_reg$LL## [1] 115.5387our_reg$vcov## (Intercept) salary job_alt difficulty
## (Intercept) 1.463992445 -0.043196792 -0.022843551 -0.109546213
## salary -0.043196792 0.004443697 0.004609611 -0.006737603
## job_alt -0.022843551 0.004609611 0.036081397 -0.024567132
## difficulty -0.109546213 -0.006737603 -0.024567132 0.073927355
## flexibility -0.011894576 0.002139895 0.003127738 -0.009827171
## market_sal 0.001509651 -0.002881612 -0.004672453 0.007909325
## flexibility market_sal
## (Intercept) -0.011894576 0.001509651
## salary 0.002139895 -0.002881612
## job_alt 0.003127738 -0.004672453
## difficulty -0.009827171 0.007909325
## flexibility 0.018456914 -0.002448031
## market_sal -0.002448031 0.002751144How Does Our Model Compare with glm?
Lets fit a model using the standard glm function and compare it to our results. Below, I output the coefficients and standard errors from both model. Not bad, right?
their_reg<-glm(decision~., data = dec%>%select(salary:decision), family = "binomial")
data.frame(our_coefs = as.vector(our_reg$betas),
our_se = sqrt(diag(our_reg$vcov)),
their_coefs = coef(their_reg),
their_se = sqrt(diag(vcov(their_reg))))%>%
kable(format = "html")%>%
kableExtra::kable_styling(bootstrap_options = "striped")| our_coefs | our_se | their_coefs | their_se | |
|---|---|---|---|---|
| (Intercept) | -0.1935129 | 1.2099556 | -0.1935129 | 1.2099512 |
| salary | 0.3384038 | 0.0666611 | 0.3384038 | 0.0666605 |
| job_alt | 1.0078720 | 0.1899510 | 1.0078720 | 0.1899494 |
| difficulty | -1.8396173 | 0.2718959 | -1.8396174 | 0.2718920 |
| flexibility | 0.2497007 | 0.1358562 | 0.2497007 | 0.1358555 |
| market_sal | -0.2805268 | 0.0524514 | -0.2805268 | 0.0524508 |
Bayesian Approach
Our maximum likelihood implementation is nice, but it has some important drawbacks. First it does not take into account prior information making it difficult to update our model as more information comes is. Second, it fails to acknowledge the uncertainty in our parameter estimates in a meaningful way. The bayesian perspective addresses both of these short commings by encorporating prior information into our understanding of observed events. The bayesian perspective is premised on bayes rule which states
\[p(A|B) = \frac{p(B|A)p(A)}{p(B)}\]
Let’s put things in terms of our logistic regression model. We have a set of parameters, \(\beta\), our data, \(D\), and a set of prior beliefs about beta \(\beta_{prior}\). We are interested in describing the density of \(p(\beta|D)\). However, the maximum likelihood approach we solved for aboove maximizes focuses exclusively on \(p(D|\beta)\). We can implement bayes rule in this scenario by modeling
\[p(\beta|D) = \frac{p(D|\beta)p(\beta|\beta_{prior})}{p(D)}\]
Lucky for us there are some familiar terms! We know that the likelihood formula above is equivalent to the \(p(D|\beta)\). The term \(p(\beta|\beta_{prior})\) is simply a probability density evaluated at \(\beta\) given the prior parameters. For example, if we use a noral prior on the regression coefficients, \(p(\beta|\beta_{prior})\) could be solved for using dmvnorm($\beta$,.
\[p(\beta|D) = \frac{L(\beta)p(\beta|\beta_{prior})}{p(D)}\]
Where \(L(\beta)\) is the likelihood of the observed data given beta, \(p(\beta)\) is equal to the probability of beta given our prior beliefs, and \(p(D)\) is a scaling constant that ensures the product of the probability probability density integrates to 0.
How can we solve for the probability of beta given our prior beliefs?
Bayesian regression has been written about extensively elsewhere. For two great articles see the blog post by count bayesie here and a more applied perspective here.
library(brms)## Loading required package: Rcpp## Registered S3 method overwritten by 'xts':
## method from
## as.zoo.xts zoo## Loading 'brms' package (version 2.10.0). Useful instructions
## can be found by typing help('brms'). A more detailed introduction
## to the package is available through vignette('brms_overview').m1priors <- c(
prior(normal(0,3), class = "Intercept"),
prior(normal(0,3), class = "b", coef = "flexibility"),
prior(normal(0,3), class = "b", coef = "job_alt"),
prior(normal(0,3), class = "b", coef = "salary"),
prior(normal(0,3), class = "b", coef = "difficulty"),
prior(normal(0,3), class = "b", coef = "market_sal")
)
m1 <- brms::brm(
decision ~ salary+job_alt+difficulty+flexibility+market_sal,
data = dec,
prior = m1priors,
family = "bernoulli",
seed = 123 # Adding a seed makes results reproducible.
)## Compiling the C++ model## Start sampling##
## SAMPLING FOR MODEL 'ce4998d23f88831d1134f9aaae974e44' NOW (CHAIN 1).
## Chain 1:
## Chain 1: Gradient evaluation took 0 seconds
## Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0 seconds.
## Chain 1: Adjust your expectations accordingly!
## Chain 1:
## Chain 1:
## Chain 1: Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 1: Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 1: Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 1: Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 1: Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 1: Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 1: Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 1: Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 1: Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 1: Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 1: Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 1: Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 1:
## Chain 1: Elapsed Time: 0.15 seconds (Warm-up)
## Chain 1: 0.151 seconds (Sampling)
## Chain 1: 0.301 seconds (Total)
## Chain 1:
##
## SAMPLING FOR MODEL 'ce4998d23f88831d1134f9aaae974e44' NOW (CHAIN 2).
## Chain 2:
## Chain 2: Gradient evaluation took 0 seconds
## Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 0 seconds.
## Chain 2: Adjust your expectations accordingly!
## Chain 2:
## Chain 2:
## Chain 2: Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 2: Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 2: Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 2: Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 2: Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 2: Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 2: Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 2: Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 2: Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 2: Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 2: Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 2: Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 2:
## Chain 2: Elapsed Time: 0.155 seconds (Warm-up)
## Chain 2: 0.137 seconds (Sampling)
## Chain 2: 0.292 seconds (Total)
## Chain 2:
##
## SAMPLING FOR MODEL 'ce4998d23f88831d1134f9aaae974e44' NOW (CHAIN 3).
## Chain 3:
## Chain 3: Gradient evaluation took 0 seconds
## Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 0 seconds.
## Chain 3: Adjust your expectations accordingly!
## Chain 3:
## Chain 3:
## Chain 3: Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 3: Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 3: Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 3: Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 3: Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 3: Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 3: Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 3: Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 3: Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 3: Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 3: Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 3: Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 3:
## Chain 3: Elapsed Time: 0.151 seconds (Warm-up)
## Chain 3: 0.154 seconds (Sampling)
## Chain 3: 0.305 seconds (Total)
## Chain 3:
##
## SAMPLING FOR MODEL 'ce4998d23f88831d1134f9aaae974e44' NOW (CHAIN 4).
## Chain 4:
## Chain 4: Gradient evaluation took 0 seconds
## Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 0 seconds.
## Chain 4: Adjust your expectations accordingly!
## Chain 4:
## Chain 4:
## Chain 4: Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 4: Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 4: Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 4: Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 4: Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 4: Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 4: Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 4: Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 4: Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 4: Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 4: Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 4: Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 4:
## Chain 4: Elapsed Time: 0.176 seconds (Warm-up)
## Chain 4: 0.152 seconds (Sampling)
## Chain 4: 0.328 seconds (Total)
## Chain 4:summary(m1)## Family: bernoulli
## Links: mu = logit
## Formula: decision ~ salary + job_alt + difficulty + flexibility + market_sal
## Data: dec (Number of observations: 250)
## Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
## total post-warmup samples = 4000
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept -0.27 1.24 -2.70 2.08 1.00 4339 3269
## salary 0.35 0.07 0.22 0.49 1.00 2138 2672
## job_alt 1.03 0.20 0.65 1.43 1.00 2610 3009
## difficulty -1.89 0.28 -2.46 -1.39 1.00 2348 2873
## flexibility 0.26 0.14 -0.00 0.54 1.00 2852 2742
## market_sal -0.29 0.05 -0.40 -0.19 1.00 1851 2607
##
## Samples were drawn using sampling(NUTS). For each parameter, Eff.Sample
## is a crude measure of effective sample size, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).vcov(m1)## Intercept salary job_alt difficulty
## Intercept 1.531707288 -0.046740837 -0.023301793 -0.110824833
## salary -0.046740837 0.004661783 0.005029606 -0.006884337
## job_alt -0.023301793 0.005029606 0.038938794 -0.027033061
## difficulty -0.110824833 -0.006884337 -0.027033061 0.076852588
## flexibility -0.013415520 0.002332528 0.002217166 -0.009911233
## market_sal 0.002752573 -0.002995312 -0.004999344 0.008067298
## flexibility market_sal
## Intercept -0.013415520 0.002752573
## salary 0.002332528 -0.002995312
## job_alt 0.002217166 -0.004999344
## difficulty -0.009911233 0.008067298
## flexibility 0.019457293 -0.002571221
## market_sal -0.002571221 0.002831588library(tidybayes)
get_variables(m1)## [1] "b_Intercept" "b_salary" "b_job_alt" "b_difficulty"
## [5] "b_flexibility" "b_market_sal" "lp__" "accept_stat__"
## [9] "stepsize__" "treedepth__" "n_leapfrog__" "divergent__"
## [13] "energy__"post<-m1%>%spread_draws(b_Intercept, b_salary, b_job_alt, b_difficulty, b_flexibility, b_market_sal)
post<-post%>%
select(b_Intercept:b_market_sal)%>%
gather(key = var, value = value)%>%
group_by(var)%>%
summarise(mean = mean(value),
sd = sd(value))fe<-c("Intercept" = 0, "cost" = .3, "job_alt" = 1, "flexibility" = .3, "difficulty" = -2, "market_sal" = -.15, sal_by_market_sal = -.003)
iter<-10
# Initialize List
model_list<-list()
model_list[[1]]<-post
for (j in 2:iter){
# Simulate People --------------------------------------------------------------------------------
# Simulate People
people<-suppressMessages(sim_people(fixed_effects = fe, n_people = 200, n_blocks = 50, p_blocks = rep(.02, 50)))
#########################################
######## Simulating decisions ###########
#########################################
# The following code chunch has the simulated people (i.e., job applicants) choose between accepting and rejecting the simulated job offers
# I then organize the data frame a bit more so that the meaning of variables is clear
dec<-suppressMessages(sim_dec(design = design_df, people = people))
dec<-dec%>%
ungroup()%>%
filter(alt == 1)%>%
select(job_id = BLOCK, applicant_id = id, salary = cost, job_alt:market_sal, decision)
# Set Priors -----------------------------------------------------------------------------------
m2prior_df <- post%>%
transmute(prior = paste0("normal(", mean, ",", sd, ")"),
class = if_else(str_detect(var, "Intercept"), "Intercept", "b"),
coef = if_else(str_detect(var, "Intercept"), "", str_remove(var, "b_")))
beta_df<-m2prior_df%>%
filter(class != "Intercept")
intercept_df<-m2prior_df%>%
filter(class == "Intercept")
m2_priors<-empty_prior()
for(i in 1:nrow(beta_df)){
m2_priors[i,]<- prior_string(beta_df$prior[i], class = beta_df$class[i], coef = beta_df$coef[i])
}
m2_priors[nrow(m2prior_df),]<- prior_string(intercept_df$prior[1], class = intercept_df$class[1])
m2 <- brms::brm(
decision ~ salary+job_alt+difficulty+flexibility+market_sal,
data = dec,
prior = m2_priors,
family = "bernoulli",
seed = 123 # Adding a seed makes results reproducible.
)
post<-m2%>%spread_draws(b_Intercept, b_salary, b_job_alt, b_difficulty, b_flexibility, b_market_sal)
post<-post%>%
select(b_Intercept:b_market_sal)%>%
gather(key = var, value = value)%>%
group_by(var)%>%
summarise(mean = mean(value),
sd = sd(value))%>%
mutate(iteration = j)
model_list[[j]]<-post
}## Compiling the C++ model## Start sampling##
## SAMPLING FOR MODEL 'e5914f8bfac5321ebcb0463913bbd85c' NOW (CHAIN 1).
## Chain 1:
## Chain 1: Gradient evaluation took 0 seconds
## Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0 seconds.
## Chain 1: Adjust your expectations accordingly!
## Chain 1:
## Chain 1:
## Chain 1: Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 1: Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 1: Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 1: Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 1: Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 1: Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 1: Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 1: Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 1: Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 1: Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 1: Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 1: Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 1:
## Chain 1: Elapsed Time: 0.103 seconds (Warm-up)
## Chain 1: 0.074 seconds (Sampling)
## Chain 1: 0.177 seconds (Total)
## Chain 1:
##
## SAMPLING FOR MODEL 'e5914f8bfac5321ebcb0463913bbd85c' NOW (CHAIN 2).
## Chain 2:
## Chain 2: Gradient evaluation took 0 seconds
## Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 0 seconds.
## Chain 2: Adjust your expectations accordingly!
## Chain 2:
## Chain 2:
## Chain 2: Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 2: Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 2: Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 2: Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 2: Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 2: Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 2: Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 2: Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 2: Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 2: Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 2: Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 2: Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 2:
## Chain 2: Elapsed Time: 0.1 seconds (Warm-up)
## Chain 2: 0.075 seconds (Sampling)
## Chain 2: 0.175 seconds (Total)
## Chain 2:
##
## SAMPLING FOR MODEL 'e5914f8bfac5321ebcb0463913bbd85c' NOW (CHAIN 3).
## Chain 3:
## Chain 3: Gradient evaluation took 0 seconds
## Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 0 seconds.
## Chain 3: Adjust your expectations accordingly!
## Chain 3:
## Chain 3:
## Chain 3: Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 3: Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 3: Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 3: Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 3: Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 3: Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 3: Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 3: Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 3: Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 3: Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 3: Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 3: Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 3:
## Chain 3: Elapsed Time: 0.099 seconds (Warm-up)
## Chain 3: 0.078 seconds (Sampling)
## Chain 3: 0.177 seconds (Total)
## Chain 3:
##
## SAMPLING FOR MODEL 'e5914f8bfac5321ebcb0463913bbd85c' NOW (CHAIN 4).
## Chain 4:
## Chain 4: Gradient evaluation took 0 seconds
## Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 0 seconds.
## Chain 4: Adjust your expectations accordingly!
## Chain 4:
## Chain 4:
## Chain 4: Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 4: Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 4: Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 4: Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 4: Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 4: Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 4: Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 4: Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 4: Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 4: Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 4: Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 4: Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 4:
## Chain 4: Elapsed Time: 0.097 seconds (Warm-up)
## Chain 4: 0.072 seconds (Sampling)
## Chain 4: 0.169 seconds (Total)
## Chain 4:## Compiling the C++ model
## Start sampling##
## SAMPLING FOR MODEL '5acd3e8e70332bf934c07cbdede94b35' NOW (CHAIN 1).
## Chain 1:
## Chain 1: Gradient evaluation took 0.001 seconds
## Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 10 seconds.
## Chain 1: Adjust your expectations accordingly!
## Chain 1:
## Chain 1:
## Chain 1: Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 1: Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 1: Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 1: Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 1: Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 1: Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 1: Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 1: Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 1: Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 1: Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 1: Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 1: Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 1:
## Chain 1: Elapsed Time: 0.093 seconds (Warm-up)
## Chain 1: 0.073 seconds (Sampling)
## Chain 1: 0.166 seconds (Total)
## Chain 1:
##
## SAMPLING FOR MODEL '5acd3e8e70332bf934c07cbdede94b35' NOW (CHAIN 2).
## Chain 2:
## Chain 2: Gradient evaluation took 0 seconds
## Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 0 seconds.
## Chain 2: Adjust your expectations accordingly!
## Chain 2:
## Chain 2:
## Chain 2: Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 2: Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 2: Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 2: Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 2: Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 2: Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 2: Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 2: Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 2: Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 2: Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 2: Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 2: Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 2:
## Chain 2: Elapsed Time: 0.09 seconds (Warm-up)
## Chain 2: 0.074 seconds (Sampling)
## Chain 2: 0.164 seconds (Total)
## Chain 2:
##
## SAMPLING FOR MODEL '5acd3e8e70332bf934c07cbdede94b35' NOW (CHAIN 3).
## Chain 3:
## Chain 3: Gradient evaluation took 0 seconds
## Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 0 seconds.
## Chain 3: Adjust your expectations accordingly!
## Chain 3:
## Chain 3:
## Chain 3: Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 3: Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 3: Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 3: Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 3: Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 3: Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 3: Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 3: Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 3: Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 3: Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 3: Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 3: Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 3:
## Chain 3: Elapsed Time: 0.09 seconds (Warm-up)
## Chain 3: 0.074 seconds (Sampling)
## Chain 3: 0.164 seconds (Total)
## Chain 3:
##
## SAMPLING FOR MODEL '5acd3e8e70332bf934c07cbdede94b35' NOW (CHAIN 4).
## Chain 4:
## Chain 4: Gradient evaluation took 0 seconds
## Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 0 seconds.
## Chain 4: Adjust your expectations accordingly!
## Chain 4:
## Chain 4:
## Chain 4: Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 4: Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 4: Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 4: Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 4: Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 4: Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 4: Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 4: Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 4: Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 4: Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 4: Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 4: Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 4:
## Chain 4: Elapsed Time: 0.092 seconds (Warm-up)
## Chain 4: 0.072 seconds (Sampling)
## Chain 4: 0.164 seconds (Total)
## Chain 4:## Compiling the C++ model
## Start sampling##
## SAMPLING FOR MODEL 'f86eff7d8bccfa331345c4eb3e100f77' NOW (CHAIN 1).
## Chain 1:
## Chain 1: Gradient evaluation took 0 seconds
## Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0 seconds.
## Chain 1: Adjust your expectations accordingly!
## Chain 1:
## Chain 1:
## Chain 1: Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 1: Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 1: Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 1: Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 1: Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 1: Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 1: Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 1: Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 1: Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 1: Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 1: Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 1: Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 1:
## Chain 1: Elapsed Time: 0.089 seconds (Warm-up)
## Chain 1: 0.07 seconds (Sampling)
## Chain 1: 0.159 seconds (Total)
## Chain 1:
##
## SAMPLING FOR MODEL 'f86eff7d8bccfa331345c4eb3e100f77' NOW (CHAIN 2).
## Chain 2:
## Chain 2: Gradient evaluation took 0 seconds
## Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 0 seconds.
## Chain 2: Adjust your expectations accordingly!
## Chain 2:
## Chain 2:
## Chain 2: Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 2: Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 2: Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 2: Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 2: Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 2: Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 2: Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 2: Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 2: Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 2: Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 2: Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 2: Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 2:
## Chain 2: Elapsed Time: 0.085 seconds (Warm-up)
## Chain 2: 0.084 seconds (Sampling)
## Chain 2: 0.169 seconds (Total)
## Chain 2:
##
## SAMPLING FOR MODEL 'f86eff7d8bccfa331345c4eb3e100f77' NOW (CHAIN 3).
## Chain 3:
## Chain 3: Gradient evaluation took 0 seconds
## Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 0 seconds.
## Chain 3: Adjust your expectations accordingly!
## Chain 3:
## Chain 3:
## Chain 3: Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 3: Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 3: Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 3: Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 3: Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 3: Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 3: Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 3: Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 3: Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 3: Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 3: Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 3: Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 3:
## Chain 3: Elapsed Time: 0.085 seconds (Warm-up)
## Chain 3: 0.095 seconds (Sampling)
## Chain 3: 0.18 seconds (Total)
## Chain 3:
##
## SAMPLING FOR MODEL 'f86eff7d8bccfa331345c4eb3e100f77' NOW (CHAIN 4).
## Chain 4:
## Chain 4: Gradient evaluation took 0 seconds
## Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 0 seconds.
## Chain 4: Adjust your expectations accordingly!
## Chain 4:
## Chain 4:
## Chain 4: Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 4: Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 4: Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 4: Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 4: Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 4: Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 4: Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 4: Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 4: Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 4: Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 4: Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 4: Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 4:
## Chain 4: Elapsed Time: 0.092 seconds (Warm-up)
## Chain 4: 0.07 seconds (Sampling)
## Chain 4: 0.162 seconds (Total)
## Chain 4:## Compiling the C++ model
## Start sampling##
## SAMPLING FOR MODEL 'e6e65c96a14556ab24bd0019d53f6429' NOW (CHAIN 1).
## Chain 1:
## Chain 1: Gradient evaluation took 0.001 seconds
## Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 10 seconds.
## Chain 1: Adjust your expectations accordingly!
## Chain 1:
## Chain 1:
## Chain 1: Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 1: Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 1: Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 1: Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 1: Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 1: Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 1: Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 1: Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 1: Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 1: Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 1: Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 1: Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 1:
## Chain 1: Elapsed Time: 0.088 seconds (Warm-up)
## Chain 1: 0.08 seconds (Sampling)
## Chain 1: 0.168 seconds (Total)
## Chain 1:
##
## SAMPLING FOR MODEL 'e6e65c96a14556ab24bd0019d53f6429' NOW (CHAIN 2).
## Chain 2:
## Chain 2: Gradient evaluation took 0 seconds
## Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 0 seconds.
## Chain 2: Adjust your expectations accordingly!
## Chain 2:
## Chain 2:
## Chain 2: Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 2: Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 2: Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 2: Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 2: Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 2: Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 2: Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 2: Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 2: Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 2: Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 2: Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 2: Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 2:
## Chain 2: Elapsed Time: 0.088 seconds (Warm-up)
## Chain 2: 0.075 seconds (Sampling)
## Chain 2: 0.163 seconds (Total)
## Chain 2:
##
## SAMPLING FOR MODEL 'e6e65c96a14556ab24bd0019d53f6429' NOW (CHAIN 3).
## Chain 3:
## Chain 3: Gradient evaluation took 0 seconds
## Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 0 seconds.
## Chain 3: Adjust your expectations accordingly!
## Chain 3:
## Chain 3:
## Chain 3: Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 3: Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 3: Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 3: Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 3: Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 3: Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 3: Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 3: Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 3: Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 3: Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 3: Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 3: Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 3:
## Chain 3: Elapsed Time: 0.091 seconds (Warm-up)
## Chain 3: 0.087 seconds (Sampling)
## Chain 3: 0.178 seconds (Total)
## Chain 3:
##
## SAMPLING FOR MODEL 'e6e65c96a14556ab24bd0019d53f6429' NOW (CHAIN 4).
## Chain 4:
## Chain 4: Gradient evaluation took 0 seconds
## Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 0 seconds.
## Chain 4: Adjust your expectations accordingly!
## Chain 4:
## Chain 4:
## Chain 4: Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 4: Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 4: Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 4: Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 4: Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 4: Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 4: Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 4: Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 4: Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 4: Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 4: Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 4: Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 4:
## Chain 4: Elapsed Time: 0.087 seconds (Warm-up)
## Chain 4: 0.071 seconds (Sampling)
## Chain 4: 0.158 seconds (Total)
## Chain 4:## Compiling the C++ model
## Start sampling##
## SAMPLING FOR MODEL '5e74c54776c4b0bd7027235790cc2276' NOW (CHAIN 1).
## Chain 1:
## Chain 1: Gradient evaluation took 0 seconds
## Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0 seconds.
## Chain 1: Adjust your expectations accordingly!
## Chain 1:
## Chain 1:
## Chain 1: Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 1: Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 1: Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 1: Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 1: Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 1: Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 1: Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 1: Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 1: Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 1: Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 1: Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 1: Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 1:
## Chain 1: Elapsed Time: 0.092 seconds (Warm-up)
## Chain 1: 0.106 seconds (Sampling)
## Chain 1: 0.198 seconds (Total)
## Chain 1:
##
## SAMPLING FOR MODEL '5e74c54776c4b0bd7027235790cc2276' NOW (CHAIN 2).
## Chain 2:
## Chain 2: Gradient evaluation took 0 seconds
## Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 0 seconds.
## Chain 2: Adjust your expectations accordingly!
## Chain 2:
## Chain 2:
## Chain 2: Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 2: Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 2: Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 2: Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 2: Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 2: Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 2: Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 2: Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 2: Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 2: Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 2: Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 2: Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 2:
## Chain 2: Elapsed Time: 0.089 seconds (Warm-up)
## Chain 2: 0.11 seconds (Sampling)
## Chain 2: 0.199 seconds (Total)
## Chain 2:
##
## SAMPLING FOR MODEL '5e74c54776c4b0bd7027235790cc2276' NOW (CHAIN 3).
## Chain 3:
## Chain 3: Gradient evaluation took 0 seconds
## Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 0 seconds.
## Chain 3: Adjust your expectations accordingly!
## Chain 3:
## Chain 3:
## Chain 3: Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 3: Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 3: Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 3: Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 3: Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 3: Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 3: Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 3: Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 3: Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 3: Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 3: Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 3: Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 3:
## Chain 3: Elapsed Time: 0.088 seconds (Warm-up)
## Chain 3: 0.071 seconds (Sampling)
## Chain 3: 0.159 seconds (Total)
## Chain 3:
##
## SAMPLING FOR MODEL '5e74c54776c4b0bd7027235790cc2276' NOW (CHAIN 4).
## Chain 4:
## Chain 4: Gradient evaluation took 0 seconds
## Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 0 seconds.
## Chain 4: Adjust your expectations accordingly!
## Chain 4:
## Chain 4:
## Chain 4: Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 4: Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 4: Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 4: Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 4: Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 4: Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 4: Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 4: Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 4: Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 4: Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 4: Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 4: Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 4:
## Chain 4: Elapsed Time: 0.087 seconds (Warm-up)
## Chain 4: 0.07 seconds (Sampling)
## Chain 4: 0.157 seconds (Total)
## Chain 4:## Compiling the C++ model
## Start sampling##
## SAMPLING FOR MODEL 'd21f4c7942fc20f9d7448bae96e30acb' NOW (CHAIN 1).
## Chain 1:
## Chain 1: Gradient evaluation took 0 seconds
## Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0 seconds.
## Chain 1: Adjust your expectations accordingly!
## Chain 1:
## Chain 1:
## Chain 1: Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 1: Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 1: Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 1: Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 1: Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 1: Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 1: Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 1: Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 1: Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 1: Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 1: Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 1: Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 1:
## Chain 1: Elapsed Time: 0.086 seconds (Warm-up)
## Chain 1: 0.102 seconds (Sampling)
## Chain 1: 0.188 seconds (Total)
## Chain 1:
##
## SAMPLING FOR MODEL 'd21f4c7942fc20f9d7448bae96e30acb' NOW (CHAIN 2).
## Chain 2:
## Chain 2: Gradient evaluation took 0 seconds
## Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 0 seconds.
## Chain 2: Adjust your expectations accordingly!
## Chain 2:
## Chain 2:
## Chain 2: Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 2: Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 2: Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 2: Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 2: Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 2: Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 2: Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 2: Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 2: Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 2: Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 2: Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 2: Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 2:
## Chain 2: Elapsed Time: 0.089 seconds (Warm-up)
## Chain 2: 0.094 seconds (Sampling)
## Chain 2: 0.183 seconds (Total)
## Chain 2:
##
## SAMPLING FOR MODEL 'd21f4c7942fc20f9d7448bae96e30acb' NOW (CHAIN 3).
## Chain 3:
## Chain 3: Gradient evaluation took 0 seconds
## Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 0 seconds.
## Chain 3: Adjust your expectations accordingly!
## Chain 3:
## Chain 3:
## Chain 3: Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 3: Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 3: Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 3: Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 3: Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 3: Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 3: Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 3: Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 3: Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 3: Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 3: Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 3: Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 3:
## Chain 3: Elapsed Time: 0.083 seconds (Warm-up)
## Chain 3: 0.069 seconds (Sampling)
## Chain 3: 0.152 seconds (Total)
## Chain 3:
##
## SAMPLING FOR MODEL 'd21f4c7942fc20f9d7448bae96e30acb' NOW (CHAIN 4).
## Chain 4:
## Chain 4: Gradient evaluation took 0 seconds
## Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 0 seconds.
## Chain 4: Adjust your expectations accordingly!
## Chain 4:
## Chain 4:
## Chain 4: Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 4: Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 4: Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 4: Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 4: Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 4: Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 4: Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 4: Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 4: Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 4: Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 4: Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 4: Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 4:
## Chain 4: Elapsed Time: 0.092 seconds (Warm-up)
## Chain 4: 0.119 seconds (Sampling)
## Chain 4: 0.211 seconds (Total)
## Chain 4:## Compiling the C++ model
## Start sampling##
## SAMPLING FOR MODEL 'd7da2b47ee98e0648c3857f75c3cfae8' NOW (CHAIN 1).
## Chain 1:
## Chain 1: Gradient evaluation took 0 seconds
## Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0 seconds.
## Chain 1: Adjust your expectations accordingly!
## Chain 1:
## Chain 1:
## Chain 1: Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 1: Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 1: Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 1: Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 1: Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 1: Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 1: Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 1: Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 1: Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 1: Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 1: Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 1: Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 1:
## Chain 1: Elapsed Time: 0.091 seconds (Warm-up)
## Chain 1: 0.104 seconds (Sampling)
## Chain 1: 0.195 seconds (Total)
## Chain 1:
##
## SAMPLING FOR MODEL 'd7da2b47ee98e0648c3857f75c3cfae8' NOW (CHAIN 2).
## Chain 2:
## Chain 2: Gradient evaluation took 0 seconds
## Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 0 seconds.
## Chain 2: Adjust your expectations accordingly!
## Chain 2:
## Chain 2:
## Chain 2: Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 2: Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 2: Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 2: Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 2: Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 2: Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 2: Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 2: Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 2: Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 2: Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 2: Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 2: Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 2:
## Chain 2: Elapsed Time: 0.092 seconds (Warm-up)
## Chain 2: 0.088 seconds (Sampling)
## Chain 2: 0.18 seconds (Total)
## Chain 2:
##
## SAMPLING FOR MODEL 'd7da2b47ee98e0648c3857f75c3cfae8' NOW (CHAIN 3).
## Chain 3:
## Chain 3: Gradient evaluation took 0 seconds
## Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 0 seconds.
## Chain 3: Adjust your expectations accordingly!
## Chain 3:
## Chain 3:
## Chain 3: Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 3: Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 3: Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 3: Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 3: Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 3: Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 3: Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 3: Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 3: Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 3: Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 3: Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 3: Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 3:
## Chain 3: Elapsed Time: 0.09 seconds (Warm-up)
## Chain 3: 0.067 seconds (Sampling)
## Chain 3: 0.157 seconds (Total)
## Chain 3:
##
## SAMPLING FOR MODEL 'd7da2b47ee98e0648c3857f75c3cfae8' NOW (CHAIN 4).
## Chain 4:
## Chain 4: Gradient evaluation took 0 seconds
## Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 0 seconds.
## Chain 4: Adjust your expectations accordingly!
## Chain 4:
## Chain 4:
## Chain 4: Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 4: Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 4: Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 4: Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 4: Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 4: Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 4: Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 4: Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 4: Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 4: Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 4: Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 4: Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 4:
## Chain 4: Elapsed Time: 0.092 seconds (Warm-up)
## Chain 4: 0.088 seconds (Sampling)
## Chain 4: 0.18 seconds (Total)
## Chain 4:## Compiling the C++ model
## Start sampling##
## SAMPLING FOR MODEL '32eb1f15a7e4a2fb557d9b4902c9fc22' NOW (CHAIN 1).
## Chain 1:
## Chain 1: Gradient evaluation took 0 seconds
## Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0 seconds.
## Chain 1: Adjust your expectations accordingly!
## Chain 1:
## Chain 1:
## Chain 1: Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 1: Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 1: Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 1: Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 1: Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 1: Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 1: Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 1: Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 1: Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 1: Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 1: Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 1: Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 1:
## Chain 1: Elapsed Time: 0.092 seconds (Warm-up)
## Chain 1: 0.069 seconds (Sampling)
## Chain 1: 0.161 seconds (Total)
## Chain 1:
##
## SAMPLING FOR MODEL '32eb1f15a7e4a2fb557d9b4902c9fc22' NOW (CHAIN 2).
## Chain 2:
## Chain 2: Gradient evaluation took 0 seconds
## Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 0 seconds.
## Chain 2: Adjust your expectations accordingly!
## Chain 2:
## Chain 2:
## Chain 2: Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 2: Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 2: Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 2: Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 2: Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 2: Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 2: Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 2: Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 2: Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 2: Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 2: Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 2: Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 2:
## Chain 2: Elapsed Time: 0.091 seconds (Warm-up)
## Chain 2: 0.072 seconds (Sampling)
## Chain 2: 0.163 seconds (Total)
## Chain 2:
##
## SAMPLING FOR MODEL '32eb1f15a7e4a2fb557d9b4902c9fc22' NOW (CHAIN 3).
## Chain 3:
## Chain 3: Gradient evaluation took 0 seconds
## Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 0 seconds.
## Chain 3: Adjust your expectations accordingly!
## Chain 3:
## Chain 3:
## Chain 3: Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 3: Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 3: Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 3: Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 3: Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 3: Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 3: Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 3: Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 3: Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 3: Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 3: Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 3: Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 3:
## Chain 3: Elapsed Time: 0.089 seconds (Warm-up)
## Chain 3: 0.106 seconds (Sampling)
## Chain 3: 0.195 seconds (Total)
## Chain 3:
##
## SAMPLING FOR MODEL '32eb1f15a7e4a2fb557d9b4902c9fc22' NOW (CHAIN 4).
## Chain 4:
## Chain 4: Gradient evaluation took 0 seconds
## Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 0 seconds.
## Chain 4: Adjust your expectations accordingly!
## Chain 4:
## Chain 4:
## Chain 4: Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 4: Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 4: Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 4: Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 4: Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 4: Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 4: Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 4: Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 4: Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 4: Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 4: Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 4: Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 4:
## Chain 4: Elapsed Time: 0.093 seconds (Warm-up)
## Chain 4: 0.075 seconds (Sampling)
## Chain 4: 0.168 seconds (Total)
## Chain 4:## Compiling the C++ model
## Start sampling##
## SAMPLING FOR MODEL '86dce6d7066dea3910a5543979034dac' NOW (CHAIN 1).
## Chain 1:
## Chain 1: Gradient evaluation took 0 seconds
## Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0 seconds.
## Chain 1: Adjust your expectations accordingly!
## Chain 1:
## Chain 1:
## Chain 1: Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 1: Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 1: Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 1: Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 1: Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 1: Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 1: Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 1: Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 1: Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 1: Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 1: Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 1: Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 1:
## Chain 1: Elapsed Time: 0.093 seconds (Warm-up)
## Chain 1: 0.113 seconds (Sampling)
## Chain 1: 0.206 seconds (Total)
## Chain 1:
##
## SAMPLING FOR MODEL '86dce6d7066dea3910a5543979034dac' NOW (CHAIN 2).
## Chain 2:
## Chain 2: Gradient evaluation took 0 seconds
## Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 0 seconds.
## Chain 2: Adjust your expectations accordingly!
## Chain 2:
## Chain 2:
## Chain 2: Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 2: Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 2: Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 2: Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 2: Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 2: Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 2: Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 2: Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 2: Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 2: Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 2: Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 2: Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 2:
## Chain 2: Elapsed Time: 0.09 seconds (Warm-up)
## Chain 2: 0.067 seconds (Sampling)
## Chain 2: 0.157 seconds (Total)
## Chain 2:
##
## SAMPLING FOR MODEL '86dce6d7066dea3910a5543979034dac' NOW (CHAIN 3).
## Chain 3:
## Chain 3: Gradient evaluation took 0 seconds
## Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 0 seconds.
## Chain 3: Adjust your expectations accordingly!
## Chain 3:
## Chain 3:
## Chain 3: Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 3: Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 3: Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 3: Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 3: Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 3: Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 3: Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 3: Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 3: Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 3: Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 3: Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 3: Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 3:
## Chain 3: Elapsed Time: 0.09 seconds (Warm-up)
## Chain 3: 0.121 seconds (Sampling)
## Chain 3: 0.211 seconds (Total)
## Chain 3:
##
## SAMPLING FOR MODEL '86dce6d7066dea3910a5543979034dac' NOW (CHAIN 4).
## Chain 4:
## Chain 4: Gradient evaluation took 0.001 seconds
## Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 10 seconds.
## Chain 4: Adjust your expectations accordingly!
## Chain 4:
## Chain 4:
## Chain 4: Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 4: Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 4: Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 4: Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 4: Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 4: Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 4: Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 4: Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 4: Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 4: Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 4: Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 4: Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 4:
## Chain 4: Elapsed Time: 0.092 seconds (Warm-up)
## Chain 4: 0.069 seconds (Sampling)
## Chain 4: 0.161 seconds (Total)
## Chain 4:model_df<-model_list%>%bind_rows()%>%
mutate(iteration = if_else(is.na(iteration), 1, as.double(iteration)))
ggplot(model_df, aes(x = iteration, y = mean, color = var))+
geom_line()
ggplot(model_df, aes(x = iteration, y = sd, color = var))+
geom_line()
LaPlace Approximation
# Initializing Parameters ------------------------------------------------------------------
# Define Priors and starting values
w_0<-m_0<-rep(0, ncol(X))
# Define precision
S_0<-diag(rep(5, times = ncol(X)))
# Defining Objective Function ----------------------------------------------------------------
post_ll<-function(theta){
v_hat<-v(X = X, beta = theta)
p<-p_hat(v_hat)
ll_x<-sum(y*p+(1-y)*(1-p))
ll_p<--1/2*t(theta-m_0)%*%solve(S_0)%*%(theta-m_0)
ll_post<- ll_x+ll_p
-ll_post
}
# Solving for theta --------------------------------------------------------------------------
optim<-optim(w_0, post_ll)
rs<-optim(par = w_0, post_ll)
theta<-rs$par
# Solving for Hessian -------------------------------------------------------------------------
posterior_hessian<- function(theta, X, y, sigma_0){
v_hat<-v(X, theta)
p<-p_hat(v = v_hat)
p_var<-p*(1-p)
mat_var<-list()
for(i in 1:length(p_var)){
mat_var[[i]]<- p_var[[i]]*X[i,]%*%t(X[i,])
}
ll_deriv<-Reduce("+", mat_var)
d2<-solve(sigma_0)
post_hesh <- ll_deriv+d2
solve(post_hesh)
}
p_h<-posterior_hessian(theta = theta, X = X, y = y, sigma_0 = S_0)
S_0<-p_h
m_0<-theta# Simulation Settings
fe<-c("Intercept" = 0, "cost" = .3, "job_alt" = 1, "flexibility" = .3, "difficulty" = -2, "market_sal" = -.15, sal_by_market_sal = -.003)
iter<-100
# Initialize List
coef_list<-list()
vcov_list<-list()
coef_list[[1]]<-theta
vcov_list[[1]]<-p_h
for(i in 2:iter){
# Simulate People ------------------------------------------------------------------
# Participants have their own preferences associated with each of these positons, for now we will simulate them as being fixed (i.e., an single preference parameter adequetly summarises all people)
# We simulate these preference distributions below.
people<-suppressMessages(sim_people(fixed_effects = fe, n_people = 200, n_blocks = 50, p_blocks = rep(.02, 50)))
#########################################
######## Simulating decisions ###########
#########################################
# The following code chunch has the simulated people (i.e., job applicants) choose between accepting and rejecting the simulated job offers
# I then organize the data frame a bit more so that the meaning of variables is clear
dec<-suppressMessages(sim_dec(design = design_df, people = people))
dec<-dec%>%
ungroup()%>%
filter(alt == 1)%>%
select(job_id = BLOCK, applicant_id = id, salary = cost, job_alt:market_sal, decision)
X<-dec[3:7]
y<-dec["decision"]
X<-cbind(1, X)
X<-as.matrix(X)
y<-as.matrix(y)
# Refit -------------------------------------------------------------------------------
rs<-optim(par = w_0, post_ll)
theta<-rs$par
p_h<-posterior_hessian(theta, X, y, S_0)
# Update Priors -----------------------------------------------------------------------
S_0<-p_h
m_0<-theta
# Save Models -------------------------------------------------------------------------
coef_list[[i]]<-theta
vcov_list[[i]]<-p_h
}coef_list%>%
map_dbl(., 6)%>%
plot(., type = "l")
vcov_list[[99]]## [,1] [,2] [,3] [,4]
## 1 2.778573e-02 -6.883213e-04 -5.472247e-04 -0.0022556562
## salary -6.883213e-04 6.864535e-05 8.351860e-05 -0.0001001873
## job_alt -5.472247e-04 8.351860e-05 7.213385e-04 -0.0005601686
## difficulty -2.255656e-03 -1.001873e-04 -5.601686e-04 0.0014060395
## flexibility -2.415504e-04 5.219542e-06 -3.445790e-05 -0.0000492812
## market_sal -8.024282e-05 -4.195005e-05 -6.459374e-05 0.0001153120
## [,5] [,6]
## 1 -2.415504e-04 -8.024282e-05
## salary 5.219542e-06 -4.195005e-05
## job_alt -3.445790e-05 -6.459374e-05
## difficulty -4.928120e-05 1.153120e-04
## flexibility 3.223782e-04 -1.133965e-05
## market_sal -1.133965e-05 4.108264e-05